Number-conserving theory of nuclear pairing gaps: a global assessment 
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We study odd-even mass staggering of nuclei, also called pairing gaps, using a Skyrme self- 
consistent mean-field theory and a numerically exact treatment of the pairing Hamiltonian. We find 
that the configuration-space Monte Carlo method proposed by Cerf and Martin offers a practical 
computational procedure to carry out the numerical solutions in large-dimensional model spaces. 
Refitting the global strength of the pairing interaction for 443 neutron pairing gaps in our number- 
conserving treatment, we find the correction to the pairing correlation energies and pairing gaps to 
have rms values of 0.6 MeV and 0.12 MeV, respectively. The exact treatment provides a significant 
improvement in the fit to experimental gaps, although it is partially masked by a larger rms error 
due to deficiencies in other aspects of the theory such as the mean-field energy functional. 
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I. INTRODUCTION 

Computer resources now make it possible to test the- 
ories of nuclear structure using the entire body of nu- 
clear data. One particular aspect of nuclear structure is 
pairing, which is important for determining stability and 
dynamical properties of nuclei. The Bardeen-Cooper- 
Schrieffer (BCS) theory [l| has been a paradigm for treat- 
ing nuclear pairing, but it is not well justified in finite 
nuclei. Besides its violation of particle-number conserva- 
tion, the condensate may collapse in finite systems. A 
recent global study of nuclear pairing gaps [2] found that 
~ 25% of nuclei lacked a BCS pairing condensate be- 
cause of the weakness of the interaction or a low single- 
particle density of states. The observed smoothness of 
nuclear binding energies calls for a theory that does not 
force a discontinuous jump between ground states with 
and without pairing condensates. In Ref. 2: it was found 
that a small but significant overall improvement in the- 
ory could be achieved by using the Lipkin-Nogami (LN) 
extension of BCS to correct for particle-number viola- 
tion 0. However, the LN treatment has its own limi- 
tations. For example, it becomes inaccurate near closed 
shells when implemented in the usual way 0, [H • On a 
practical level, iterative BCS-LN solvers often have con- 
vergence problems near closed shells. We note that there 
are many methods other than the LN extension of BCS 
to treat the pairing interaction more accurately, includ- 
ing direct diagonalization in truncated spaces 6]. For 
methods that emphasize particle number conservation, 
see Ref. [7[ and references therein. 

Here we address the question of the importance of a 
better treatment of pairing correlations by carrying out a 
global survey using a numerically exact technique to cal- 
culate pairing correlation energies at fixed particle num- 
ber. In particular, we employ the configuration-space 
Monte Carlo (CSMC) algorithm of Cerf and Martin Q. 
Numerically exact solutions can also be obtained by di- 
rect diagonalization of the pairing Hamiltonian in config- 



uration spaces of fixed seniority 

@,E3. but the CSMC 

method is more efficient and can be implemented in much 
larger spaces. 

Since our aim is to assess the relative performance 
of the theory with and without an exact treatment of 
pairing, we will avoid introducing extraneous elements 
and closely follow the methodology of Ref. Q. In that 
work the performance of various self-consistent mean- 
field (SCMF) methods was tested on the neutron and 
proton pairing gaps for odd- /I nuclei; here we use the 
same 443 odd-neutron gaps to assess the importance of 
an exact treatment of the pairing interaction. 

The neutron pairing gap for an (odd) neutron num- 
ber N is defined by the second-order energy difference in 
neutron number 

^ ( o\N) = -\ [E(N + 1) + E(N - 1) - 2E(N)] , (1) 

where E(N) is the ground-state energy of the nucleus 
with TV neutrons and Z protons. The proton number 
Z is the same for all three nuclei in Eq. (fTJ and is not 
indicated explicitly in the formula. 

As a prototype of an SCMF theory we use the energy 
density functional constructed from the SLy4 Skyrme 
functional for the normal density part and a density- 
dependent contact interaction for the pairing part. The 
Hartree-Fock+BCS (HF+BCS) equations are solved us- 
ing the EV8 code [12j ■ We construct a pairing Hamil- 
tonian whose single-particle orbital energies and pairing 
matrix elements are extracted from the SCMF calcula- 
tion of Ref. 0. Next, we solve this Hamiltonian exactly 
using the CSMC method, which is free of a sign prob- 
lem when all pairing matrix elements are attractive. The 
SCMF interaction energies are also taken from the cal- 
culations of Ref. 2} . The main difference here is in the 
treatment of pairing correlation energies and the strength 
of the pairing interaction. The performance of the the- 
ory is measured by the root-mean-square (rms) of the 
residuals with respect to the experimental data set after 
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making a least-squared fit of the overall pairing interac- 
tion strength. The differences between the Monte Carlo 
treatment and the BCS approximation are used to es- 
timate the importance of a particle-number-conserving 
exact treatment of pairing. 

The outline of this paper is as follows. In Sec. |TT] we 
discuss our methodology of constructing a pairing Hamil- 
tonian from the SCMF results and how we use its exact 
solution to obtain an improved estimate of the pairing 
gap. In Sec. [TIT] we describe the CSMC method used to 
solve the pairing Hamiltonian. In Sec. HVI we present our 
results for the pairing gaps. Our conclusions are given in 
Sec. El 



II. METHODOLOGY 

The leading approach in the search for a computation- 
ally tractable theory of nuclear structure starts with an 
SCMF theory to construct a set of configurations and 
then mixes the configurations through a residual interac- 
tion to restore broken symmetries and add a correlation 
contribution to the total energy. By itself, the mean-field 
theory is straightforward. However, there are different 
ways to introduce correlations, even if we limit ourselves 
to pairing correlations. The BCS treatment is the sim- 
plest way to introduce pairing correlations and can be 
easily implemented once the single-particle wave func- 
tions and energies have been obtained from the mean- 
field theory. The Hartree-Fock-Bogoliubov (HFB) ap- 
proximation is an extension that is required when the 
orbital properties depend on pairing, but like BCS it vi- 
olates particle-number conservation. To gain the benefit 
of the HFB approximation, the pairing Hamiltonian must 
be defined in very large model spaces and with a more 
general interaction than can be treated with the present 
CSMC method. Here we construct pairing Hamiltoni- 
ans in which the orbitals are fixed from the mean-field 
calculation, as in the BCS approximation, and have in- 
teraction matrix elements that are all attractive. The 
HFB might be required in the dripline region. However, 
since very few of the known experimental gaps are in this 
region, our conclusions should apply to the vast majority 
of nuclei for which data exist. 

For construction of the pairing Hamiltonian, we follow 
closely the treatment of Ref. [l3[ as implemented in EV8. 
The single-particle energies are taken directly from the 
eigenvalues of the single-particle SCMF Hamiltonian and 
the interaction is chosen as a density-dependent contact 
interaction 

V(r,r') = -V (l-V^jS(T-r') (2) 

together with an energy cutoff factor described below. 
In Eq. ^ po — 0.16 fm -3 is the conventional saturation 
density of nuclear matter and the parameter rj controls 
the specific density dependence. We will use r\ = 0.5, 



called the "mixed" density-dependent pairing interac- 
tion. The strength Vq is determined by minimizing the 
rms of the residuals of the calculated pairing gaps from 
their experimental counterparts 0]. 

As implemented in EV8, the mean field is invariant 
under time reversal and the self-consistent single-particle 
orbitals appear in degenerate time-reversed pairs i and 
i with energy e^. The total number of orbital pairs is 
Cl. The antisymmetrized pairing matrix elements Vij are 
taken to be 

ee M{a\V\jj) - (ii\V\jj))fj , (3) 

where V is given by Eq. ([2]) and /, are energy cutoff 
factors [l3| 



[i + e (e 4 -o)/6 l + e (- E «-a)/6j ^ 

with a = 5 MeV and b = 0.5 MeV. Denoting the single- 
particle orbitals by fa (r, a) , we have 

Vy = ~frfjV J dr \S2 l&M| 2 J (E l<M r V)| 2 ) 

x(l-^) -(5) 
Next we construct the pairing Hamiltonian 

n n 
H = Hi + H 2 = }^ Si(a\ai + ajaj) + 2J VyajotajOj . 
i i^j 

(6) 

Note that the Hamiltonian © does not include diagonal 
matrix elements Va. We assume that they have already 
been incorporated into the mean-field part of the energy 
density functional. 

An Hamiltonian of the form Eq. ^ was used in a re- 
cent study comparing the BCS approximation with ex- 
act matrix diagonalization results in model spaces of size 
Cl = 16 [ijj]. However the sizes of the model space re- 
quired for a global survey are prohibitively large for di- 
rect matrix diagonalization methods to be practical. We 
therefore use the CSMC method which scales much more 
gently as a function of Cl. This method can be used to 
find the exact ground-state energy -Ecsmc of the Hamil- 
tonian in Eq. ^ to within a statistical error. We note 
that the pairing Hamiltonian can be solved algebraically 
for special forms of the interaction following Richardson's 
method [l5| , but these are not applicable to more general 
interactions such as in Eq. (2). 

Our improved estimate for the total ground-state en- 
ergy is given by 

E = -Escmf — -Bbcs + Ecsmc , (7) 

where -Escmf is the SCMF energy calculated with SLy4 
plus the density-dependent contact pairing interaction 
and -Ebcs is the BCS ground-state energy of the 
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Hamiltonian in Eq. ([6]). In Eq. (J7J we are essentially re- 
placing the BCS energy of the Hamiltonian in Eq. © 
with its exact CSMC ground-state energy. Both -Escmf 
and -Ebcs are calculated with an interaction strength 
Vo determined by minimizing the rms deviation of the 
SCMF gaps from the experimental gaps. However, 
-Ecsmc is calculated with a renormalized strength, de- 
termined by minimizing the rms residuals (with respect 
to experiment) of the theoretical gaps calculated from 
Eq. Q and Eq. Q. 

One problem of using a pairing Hamiltonian from a 
theory such as the one discussed in Ref. @ is that there 
are diagonal interaction matrix elements both in the 
mean-field part as well as in the pairing part of the en- 
ergy functional. Since the pairing interaction is added to 
describe correlations beyond those obtained in the SCMF 
with a single Slater determinant, it should not add diag- 
onal interactions beyond the mean field. At an extreme, 
if the BCS condensate collapses, the BCS correlation en- 
ergy should be zero. Because of these considerations we 
do not include diagonal interaction matrix elements Vu 
in our Hamiltonian ^ for either the CSMC or the BCS 
calculations. Such matrix elements remain, however, in 
the pairing part of the SCMF theory. 



III. CONFIGURATION SPACE MONTE CARLO 
SOLVER 

The CSMC Hamiltonian solver has been applied to 
individual isotope chains [16( but our work here is its 
first use in a global survey. To introduce the various 
parameters of the method and make our presentation 
self-contained, we review the algorithm in Sec. IIII Al In 
Sec. IIII Bl we discuss the statistical Monte Carlo error 
and demonstrate the computational scaling properties of 
the method. 



paired configurations yj. Let us label the fully paired 
eigenstates of Hi, in Eq. ©, by |n) = |m, n^, . . . nn) 
where tii = or 1 is the pair occupation number of the 
i-th two-fold degenerate level, and 



ith 



ffi|n) =£ ap (n)|n> 



E sp {n) =2j2 £ i 



(9) 



(10) 



The wave function |\&(t)) can be written as a linear com- 
bination of these paired configurations |n) 



|*(r))=$> n (T)|n) 



(11) 



The coefficients a n (r) can be chosen to be all positive 
and normalized as 



1 



(12) 



Thus, the wave function |\&(t)) can be represented by an 
ensemble of paired configurations |n) that are distributed 
with probability a n (r). 

The evolution in imaginary time is carried out as a 
series of time evolutions, each of which is over a small 
time step At 



|tf (t + At)) = e -^-- Bt > AT |*(7-)) 



(13) 



Using the Suzuki- Trotter symmetric decomposition fjl\ - 
[l9j . we write the short-time propagator as 

e -(H-E t )A T _ e -(H 1 -E t )AT/2 e -H 2 Ar e -(H 1 -E t )AT/2 

+ 0{At 3 ) .(14) 



A. The Monte Carlo method 

In the following, we assume the particle number N to 
be even. The algorithm applies to Hamiltonians of the 
form in Eq. (|6]) for which all pairing interaction matrix 
elements satisfy Vij < 0. The overall operation of the 
algorithm is similar to many other Monte Carlo methods 
where a trial state |3>) is evolved in imaginary time 



|*(r)) 



-(H-B t )T| $ \ 



(8) 



where H is the system's Hamiltonian and E t is an energy 
parameter adjusted to keep the normalization of |\&(t)) 
approximately fixed. The initial evolution filters out the 
ground-state component of the initial trial state, and sub- 
sequent evolution is used to obtain better statistics for 
the ground-state energy. 

The remaining details of the Monte Carlo method relies 
on the representation of V?(t) as a superposition of pure 



We arc interested in calculating the matrix elements of 
Eq. (Tnj between two paired configurations |n) and |n'). 
Since H\ is diagonal in the |n) basis, the only non-trivial 
part is the matrix elements of e~ H2Ar . Expanding this 
propagator in a Taylor series, we have [8| 



(n'le-^^ln) 



L=0 



P{L) 



,L 



W K (n^n>) 



(15) 

where oj — ^(£1—^ + 1) and v = ujVAt defines a 
dimensionless time step with V = ^^Vij/Vt 2 . Here 

ii 

W K (n — > n') represents the weight of a path of L pair 
hops that takes the configuration n to n'. Each pair 
hop describes the transition of a pair of particles from an 
occupied two-fold level to an empty two-fold level. The 
probability to have L pair hops in the time interval At is 
a Poisson distribution P(L) — e~ u v L /LI. The parameter 
uj represents the total number of possible pair hops and 
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v is the average number of pair hops in the time interval 
At. The weight W K (n — > n') is given by 



W„(n->n') = Y[ 



m — 1 



V, 



V 



(16) 



where (i m ,i m ) and {j m ,j m ) are the orbital pairs whose 
occupations are swapped at the m-th step of the L-step 
pair hop process. 

In practice, we carry out the Monte Carlo evolution as 
follows. We take the initial state |$) to be the ground- 
state configuration of Hi and replicate it N e times to 
generate the initial ensemble. Subsequently, the mem- 
bers of the ensemble are evolved independently. For 
each time step At, the time evolution is done stochas- 
tically using Eqs. (|15l) and (TTB)) . The number of pair 
hops L is drawn from the Poisson distribution P(L), 
and an L-step pair hop process is carried out. At each 
step we choose an occupied pair orbital {i m ,im) from 
a uniform distribution, and an orbital (j m ,j m ) which 
is either unoccupied or equal to (i m ,i m )i again from 
a uniform distribution. The occupation numbers of 
(i ra ,'m) and {j m ,jm) are swapped. The resulting new 
configuration is replicated stochastically with a weight of 
cxp[-(£ sp (n) + E sp (n') - 2E t )A T /2]W K {n -> n'). 

We adjust the normalization energy E t during the time 
evolution to keep the ensemble size stable. At the fc-th 
time step we define 



E t (k) = E t (k-l) + -^-\n 
At 



iV e (fc-l) 
N e (k) 



(17) 



where N e (k) is the size of the ensemble at the fc-th time 
step. 

We repeat the above process N T times for the initial 
evolution. The resulting ensemble of fully paired 

configurations n m (m — 1, . . . , N e ) can be expressed 
as the wave function 



l*i> 



1 Nb 



(18) 



and is our first representative of the ground-state ensem- 
ble. 

The ground-state energy E can be calculated from 
the ground state \P using E = 5Z„/(n'|if|^) (where we 
have used X)n'( n 'l^ / ) = •"■)• Approximating ^ by "fi in 
Eq. ITS)) , we estimate the ground-state energy to be 



E x 



— y 

771 — 1 



[£; s p(n m ) +E v (n m )} , 



(19) 



where E sp (n) is given by Eq. (fT0|) and 



B t; (n) = ^(n'|iJ 2 |n) = ^y ij 



(20) 



The prime on the summation in Eq. (|20)) denotes that 
the sum is restricted to those combinations ij where the 
orbital pair (i, i) is occupied in |n) and the orbital pair 
is either unoccupied in |n) or the same as (i, i). 
Additional representatives of the ground-state 

wave function are generated by evolving the ensemble 
an additional number of time steps Nt and taking a 
representative every N c steps to ensure uncorrelated en- 
sembles. Using relations similar to Eq. (|T9| . we obtain 
N E = N T /N C estimators E\,E%, . . . , E Ne for the ground- 
state energy. The final estimate for the CSMC ground- 
state energy is 



1 n e 
-Ecsmc = — — > E. 

1\ F ^ — ' 



and its corresponding statistical error is 



(21) 



N E 



\ N E (N E - 1) 



Y,y,( Ei - ^CSMC) 5 



(22) 



So far, we have discussed a system with an even num- 
ber of particles N. The generalization to an odd number 
N D is straightforward. We put a single particle in one of 
the orbitals of a degenerate pair. This pair of orbitals 
becomes effectively blocked, i.e., it cannot participate in 
the pair transitions between orbitals. The energy of the 
remaining N — 1 particles is found by applying CSMC 
to the reduced space in which the blocked orbital pair is 
excluded. The total energy of the A'o-particle system is 
then given by 



E b {N ) = E b (N - 1) + e b 



(23) 



where e b is the single-particle energy of the blocked or- 
bital and E b (N — 1) is the energy of N a — 1 particle 
system in the reduced space. The ground-state energy of 
the odd-JV system is found by minimizing Eq. (|23[) over 
different choices of the blocked orbital b. 

For the calculations in this work, we have taken N T = 
5000, N T = 50, 000, N c = 500 and N e (0) = 25, 000. This 
gives N E — 100 estimators Ei of the the ground-state 
energy and its error in Eqs. (l2Tj) and (|22[) . respectively. 



B. Statistical error 

The statistical error a in the CSMC energy estimate 
can be written as 



Xe 



N E N e 



(24) 



where af n is the intrinsic variance of the energy, i.e., the 
variance of the quantity E sp (ri) + E v (n) for paired con- 
figurations n that are distributed according to a n . Here 
Ne = Nt/N c is the number of uncorrelated ensembles of 
size N e used in the CSMC calculation. The replication 
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process described in Sec. MI Al introduces correlations be- 
tween configurations in the ensemble at a given time step 
and N e /xe (with \ e > 1) represents the effective number 
of uncorrelated configurations. 

In the following, we provide an estimate for the intrin- 
sic standard deviation a m assuming a constant pairing 
interaction. In this case £^(11) in Eq. (|2"0|) is a constant, 
and af n is determined solely by the variance of the E sp (n) . 
With E sp (n) given by Eq. (fTUf , its average value E sp over 
the various configurations in the ensemble is 



where n, 



E 



n.ar, 



X)n t , (25) 



Eq. (|25j) holds for any constant A 



but we choose A to be the chemical potential to minimize 
the particle number fluctuations. 

In the Appendix we use the BCS wave function to es- 
timate the fluctuations in E sp [see Eq. (IA.4I) ]. For a uni- 



form single-particle spectrum with a bandwidth E c 
(A is the BCS pairing gap), we find 



f 



-51E r A 



» A 



(26) 



We can use this expression to estimate the scaling of a- ln 
with the size 51 of the single-particle space. For weak to 
moderate pairing, A cx 51, and <r in cx 51 3 / 2 . Our sim- 
ple estimate ([26]) is accurate to within a factor of ~ 2 
(see Fig. |5] in the Appendix). We have checked that even 
in cases when the single-particle spectrum is non- uniform 
and the pairing interaction is orbital-dependent (e.g., the 
nuclear pairing Hamiltonians used in this work), expres- 
sion ([2"fj]) (where A is taken to be an average pairing gap) 
provides a reasonable estimate for the intrinsic error in 

-Esp- 

If all iV e configurations of the ensemble at a given time 
step were to be uncorrelated, the CSMC error would have 
been a- m / V 'NeN e . Since these configurations are corre- 
lated in the CSMC calculation, the actual statistical error 
is larger by a factor of y/xe [see Eq. ([24]) ]. In Fig. Q] we 
show (solid circles) this enhancement factor (as de- 
termined empirically from the CSMC statistical error for 
v = 0.1) versus 51 for a uniform single-particle spectrum 
with level spacing of 1 MeV at half filling (N = 51) and a 
constant pairing strength of Vij = 0.3 MeV. The dashed 
line is a fit to Xe ~ 1 + (51/51o) 3 - In general we find that 
the scaling of Xe with 51 depends on the strength of the 
pairing interaction. 

To illustrate the scaling of the CSMC computational 
time with the size 51 of the single-particle space, we con- 
sider the same example as in Fig. [1] The results, shown 




FIG. 1: The factor yfxZ as a function of 51 for an equally- 
spaced single-particle spectrum with level spacing of 1 MeV 
and pairing strength of Vij = 0.3 MeV. The dashed line de- 
scribes the fit Xe = 1 + (fi/^o) 3 witn ^0 = 15.2. 



by symbols in the upper panel of Fig. [2] scale (up to an 
additive constant) as 51 2 (dashed curve). 

The lower panel of Fig. [5]shows the statistical error cal- 
culated from Eq. ([22]) . It appears to scale as 51 3 (dashed 
line) 2 . Spaces as large as 51 = 30 are easily computed and 
we shall argue below that the accuracy achieved is ade- 
quate for our purposes. In contrast, if the calculations 
were done by conventional matrix diagonalization, one 
would have to deal with a matrix of dimension 1.6 x I0 8 
and 51 = 50 (corresponding to a matrix of dimension 
1.3 x I0 14 ) would be completely out of reach. 

The CSMC calculations for the global survey were car- 
ried out using two values of v (0.05 and 0.10) and av- 
eraging their respective energy estimates. We checked 
that the corresponding time steps were sufficiently small 
to avoid a significant systematic error from the Suzuki- 
Trotter decomposition in Eq. (fill . To test for other 
biases in the CSMC algorithm, we compared with the 
matrix diagonalization of the Hamiltonian (]6]) for two 
cases presented in Ref. [14|, namely 118 Sn and 206 Pb. 
The pairing Hamiltonians were obtained from Ref. 0; 
Hi is derived from SCMF with the Skyrme SLy4 energy 
functional and H2 is of the contact form as in Eq. 
with r\ = (and no cutoff factors fi). The single-particle 
space in these examples has a size 51 = 16, which requires 
matrices of dimension ~ 13, 000 for the direct diagonal- 
ization. The calculated correlation energies (measured 
relative to the HF ground-state energy) are shown in Ta- 
ble |T] for two values of the pairing strength, Vq = 360 and 



2 This error seems to depend on the parameters of the pairing 
1 Note that fii differs from the quantum mechanical expectation Hamiltonian and for a stronger pairing interaction we find a more 

value of the pair occupation operator n;. moderate scaling of Q 2 . 



6 




FIG. 2: Scaling of computational effort with the number of 
orbitals 0, for an equally-spaced single-particle spectrum with 
level spacing of 1 MeV and a constant pairing interaction 
Vij = 0.3 MeV: a) single-processor CPU time for a single 
CSMC calculation (the dashed line is a fit describing a scal- 
ing of fi 2 ); b) statistical error of Eq. (|22[) (the dashed line 
corresponds to a scaling of fi 3 ). See text for the values of the 
CSMC parameters. 



450 MeV fm 3 . From the results we see that there are no 
discernible systematic errors in the CSMC calculations. 





Vo = 360 
Lanczos CSMC 


Vo = 450 
Lanczos CSMC 


118 Sn 

206p b 


2.564 2.569 ± 0.006 
0.363 0.365 ± 0.004 


4.553 4.546 ± 0.006 
0.626 0.626 ± 0.005 



TABLE I: Comparison of pairing correlation energies calcu- 
lated by CSMC and by exact diagonalization method (Lanc- 
zos algorithm). The interaction strength Vo is in units of 
MeV-fm 3 and the correlation energies in units of MeV. 



For the CSMC calculation in our global survey, it is im- 
portant that the Monte Carlo statistical error does not 
degrade the accuracy of the calculated pairing gaps to a 
point where the performance measure would be affected. 
The maximal permissible statistical error is estimated as 



follows. We take a typical rms of the residuals (deviations 
between theory and experiment) in the range 0.25-0.30 
MeV, and demand that the Monte Carlo statistical con- 
tribution, calculated in quadratures, be less than 0.01 
MeV. This requires that the average statistical error for 
the pairing gap be smaller than -\/0.25 2 — 0.24 2 0.07 
MeV. In fact, with our choices of the numerical param- 
eters in the global calculation, the maximal statistical 
error (~ 0.05 MeV) satisfies this upper bound for all 
cases. 



IV. RESULTS 

For our global survey of odd neutron gaps, we take the 
same nuclei as in Ref. 0, where 443 odd neutron pair- 
ing gaps were calculated and compared with experiment. 
Our procedure for obtaining a new set of theoretical gaps 
involves the following steps: 

1. We start with the full SLy4+pairing energies as cal- 
culated in Ref. 2] and construct the pairing Hamil- 
tonian H in Eq. © using the converged SCMF 
single-particle energies and wave functions. The 
diagonal interaction matrix elements Vu are not in- 
cluded in Eq. ©. 

2. We calculate the BCS ground-state energy of H, 
taking the same interaction strength Vq as in the 
original calculations (i.e., Vo = 700 MeV-fm 3 ), and 
subtract it off the total SCMF energy. 

3. We calculate the exact ground-state energy of H 
by CSMC and add it back to to obtain our new 
estimate of the ground-state energy. The overall 
interaction strength is renormalized in the CSMC 
calculation, and its value is determined by mini- 
mizing the rms residuals of the newly calculated 
pairing gaps. 

4. To make a fair comparison with the BCS, we repeat 
step 3 with BCS energy (excluding diagonal inter- 
action matrix elements as in the CSMC), refitting 
the overall strength of the interaction to minimize 
the rms residuals. 3 

Some remarks are in order regarding our refit. The 
CSMC calculations were performed at two different val- 
ues of Vq (560 and 700 MeV fm 3 ). We used a linear 
interpolation to obtain the ground-state energy for inter- 
action strengths between these two values. The new value 
of Vo is determined by minimizing the rms of the residu- 
als using a linear least-squared fit. The model space for 



3 In principle, step 4 should not be necessary but we found that 
the strength Vo reported in Ref. Q is not optimal for the BCS 
theory presented there. 
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our CSMC calculation consists of all orbitals for which 
ff > 0.01 [see Eq. (gj)]. We verified the convergence of 
our calculations by repeating them for a model space with 
ff > 0.001. The largest model space used in these calcu- 
lations is fl = 64 for the nucleus N = 156 and Z = 100. 
For this nucleus each CSMC calculation takes about 30 
minutes on a single processor. 

The results for our refits are shown in Table [TT] along 
with the fit reported in Ref. [2|. First we note that the 
fitted value of the interaction strength Vb is smaller for 
the CSMC gaps than its fitted value for the BCS gaps by 
about 6%. It is not surprising that the required strength 
is higher in a theory (e.g., BCS) that is subject to pairing 
collapse and gives a zero correlation energy in some of the 
nuclei. In fact, the differences between correlation ener- 
gies comparing the CSMC and BCS can be quite large; 
their rms difference is ~ 0.6 MeV for the more than 900 
nuclei in our data set. However, the observable quanti- 
ties are not the correlation energies but the pairing gaps, 
for which the differences are much smaller. The rms of 
the differences between the CSMC and the BCS pair- 
ing gaps is 0.12 MeV. Given that the total rms residuals 
of the theory with respect to experiment is of the order 
0.25-0.30 MeV, this difference between CSMC and BCS 
appears to be quite significant. However, one must real- 
ize that when there are independent sources of error, the 
larger ones can effectively mask the others. This can be 
seen in the third column of table II, reporting the rms 
residuals of the CSMC and BCS pairing gap with respect 
to experiment. The corresponding values, 0.28 and 0.24 
MeV, only differ by 0.04 MeV. However, this is close to 
what one would expect when adding in quadratures the 
error of the BCS approximation (0.12 MeV) and the other 
sources of error (0.24 MeV). We also note that the val- 
ues for the fitted strength and the rms residuals reported 
in Ref. Q are somewhat higher than their corresponding 
values in our BCS fit, to which it should be compared. 



Method 


Vb 


rms 




(MeV fm 3 ) 


(MeV) 


SCMF [2] 


700 


0.30 


BCS 


667 


0.28 


CSMC 


627 


0.24 



TABLE II: The rms residuals of the calculated pairing gap 
using different theoretical methods. See text for a description 
of the various methods. 



To illustrate the performance of the theory locally, we 
compare in Fig. [3] the theoretical and experimental pair- 
ing gaps for the chain of Sn isotopes (Z = 50). All three 
theories (SCMF, BCS and CSMC) overestimate the gap 
but follow correctly its overall dependence on neutron 
number, including the the dip at N — 65 and the sharp 
drop near the N — 82 shell closure. When compared 
with the experimental gaps, the CSMC shows a modest 
but systematic improvement over the SCMF and BCS 
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FIG. 3: Pairing gaps of Sn isotopes: predictions of various 
theories (SCMF, BCS and CSMC) are compared with the 
experimental gaps. 



theories, except for the N = 77 — 81 nuclei in the vicinity 
of the shell closure. 

It would be useful to know whether there are any sys- 
tematic criteria for identifying nuclei for which the im- 
proved treatment of pairing has the most benefit. One 
criterion could be the magnitude of the error (i.e., resid- 
ual) comparing the SCMF or BCS pairing gaps with their 
experimental values. To examine the dependence on the 
SCMF error we take subsets of gaps whose SCMF resid- 
ual is larger (in absolute value) than some given value 
and calculate the rms of the subset as a function of the 
lower cutoff. The results are shown in Fig. 2) The rms 
error of the subset increases with cutoff but in the CSMC 
approach it does so at a lower rate than in the BCS 
treatment. For example, when we keep only those nuclei 
whose SCMF rms residual is greater than 0.5 MeV, we 
find that the BCS rms increases to 0.68 MeV while the 
CSMC rms is only 0.52 MeV, an improvement of 0.14 
MeV. The inset of Fig. H shows the rms of the CSMC 
correction to the BCS residual versus the lower cutoff of 
the SCMF residual. This rms exhibits a gradual increase 
from about 0.12 MeV when all the nuclei are included to 
about 0.19 MeV when we include only those nuclei whose 
SCMF residual is greater than 0.5 MeV. Thus, there is a 
mild increase in the benefit derived from the exact treat- 
ment when the residual error is large. 

To narrow further the conditions under which an ex- 
act treatment is beneficial, we go back to the symmetries 
that are broken in mean field and BCS theory, namely 
particle-number conservation and rotational symmetry. 
A measure of the BCS violation of particle-number con- 
servation is given by 

(AA0 2 = {{N - {N)f = 4]T(1 - vf)vf , (27) 

where vf are the BCS occupation numbers. We divide 
the nuclei with odd number of neutrons into bins of width 
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Minimal SCMF residual of the set (MeV) 

FIG. 4: The rms residuals of the pairing gap in nuclei for 
which the absolute deviation of the SCMF pairing gap from 
the experimental value is greater than the value shown on 
the horizontal axis. The inset shows the rms of the CSMC 
correction to the BCS gap, the horizontal axis being the same 
as the main graph. 



1 according to their particle-number fluctuation AiV (nu- 
clei with AN = have their own). Fig. [5] shows the rms 
of the residuals for the nuclei in each bin versus the mid- 
point of the bin. The bin with AN = consists of all 
nuclei for which the BCS pairing has collapsed. Clearly 
the CSMC treatment is needed in that situation. The 
CSMC also gives an improvement for the bin centered at 
AiV = 3.5 having the strongest pairing condensate. This 
is likely due to the too-large pairing strength Vq required 
for the global BCS fit. Thus, when compared to the ex- 
act CSMC results, the BCS approximation seems to be 
adequate when 1 < AiV < 3. If we use in the BCS treat- 
ment the lower value of the CSMC interaction strength, 
we find the BCS performance to improve gradually with 
increasing AiV and to become comparable to the CSMC 
performance for AiV > 2. 

The violation of rotational symmetry in a nucleus is 
often characterized by the mass quadrupole deformation 
parameter fa denned by 

A" < 2s) 

where A, (r 2 ) and Qo are, respectively, the mass number, 
rms radius and intrinsic quadrupole moment of the nu- 
cleus. We divide the odd-iV nuclei into bins of width 0.1 
according to their deformation fa in the SCMF treat- 
ment. The rms of the residuals are calculated for the 
nuclei in each bin and their values (in both BCS and 
CSMC) are plotted versus the bin centers in Fig. [5] We 
observe that there is almost no difference between the two 
treatments for oblate nuclei. For spherical and strongly 
deformed prolate nuclei, the CSMC gives a moderate 
improvement over BCS while for moderately deformed 
prolate nuclei there is a significant improvement in the 
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FIG. 5: The rms residuals of the pairing gap as a function 
of particle- number fluctuation AiV [see Eq. (57)]. The nu- 
clei were divided into bins of width f according to the values 
of AiV obtained from the v amplitudes of the SCMF theory 
(AN = nuclei have their own bin). The points are posi- 
tioned at the center of the bins and the lines are drawn to 
guide the eye. 
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FIG. 6: The rms residuals of the pairing gap as a function of 
deformation /?2 [see Eq. (|28p j. The nuclei were divided into 
bins according to their /?2 value in the SCMF treatment. The 
points are positioned at the center of the bins, and the lines 
are drawn to guide the eye. 



CSMC method as compared with the BCS approxima- 
tion. 

In Fig. [7] we show the average value of the ratio 
between the CSMC correction to the BCS correlation 
energy and the CSMC correlation energy, (SEqsmg — 
SE B cs)/SE C smc, versus AiV (here SE B cs = En F -E B cs 
and SE CSMC = E HF - E CSMC ). At AiV = the BCS so- 
lution collapses to the HF solution and this ratio is just 
1. We observe the above ratio to decrease monotonically 
versus AiV as the BCS approximation becomes better. 
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FIG. 7: The ratio (SEcsmc — SEbcs) / SEcsmc as a function 
of particle-number fluctuation AiV for nuclei with odd and 
even number of neutrons N. This ratio decreases as the BCS 
approximation becomes better at larger AiV. The case AiV = 
1.5 is a borderline case: the odd-iV behaves more like AiV = 
while the even-iV is closer to AiV ^> 1. 



V. CONCLUSION 

Starting from an SCMF theory of the pairing gaps 
and treating pairing correlations exactly beyond the BCS 
approximation (with a renormalized pairing interaction 
strength) , we found a significant improvement in the the- 
ory as measured by the rms residuals of the pairing gaps. 
The exact calculations were carried out by constructing 
a pairing Hamiltonian from the SCMF output and using 
the configuration space Monte Carlo (CSMC) method. 

We find the improvement in the rms residuals of the 
pairing gaps to be most significant in nuclei for which the 
BCS condensate was weak, as measured by the smallness 
of particle-number fluctuation AiV. Based on our results, 
the BCS seems to be adequate if one limits the theory to 
nuclei for which 1 < AiV < 3. The artificially high value 
of the BCS interaction strength in a global fit leads to 
pairing gaps that are on average too large in nuclei with 
AiV > 3. We also found the improvement to be larger in 
moderately deformed (/?2 ~ 0.1 — 0.2) prolate nuclei. 

The total residual in the SCMF treatment of pairing 
gaps can be thought of as coming from two parts, the 
inadequacy in the mean field and the approximate treat- 
ment of pairing. Since in the CSMC method the pairing 
part is treated exactly, the CSMC residuals are wholly 
due to the inadequacy of the mean field. The rms of 
the residuals of the pairing gaps in the BCS treatment 
over their values in the CSMC treatment is then an es- 
timate of the error involved in the BCS approximation. 
We find this rms value to be ~ 0.12 MeV (see inset of 
Fig. 2}, and propose it as the bound on the accuracy that 
can be achieved in an SCMF theory that treats pairing 
correlations approximately 

From the computational perspective, the most notable 



aspect of this work is the use of the CSMC algorithm, 
which has not been used previously in a global survey. 
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Appendix 

In this Appendix we derive an estimate of the intrinsic 
statistical error <7i n based on the BCS wave function and 
show that it reduces to Eq. (j2"6) in the limit of a large 
bandwidth E c 3> A (assuming a uniform single-particle 
spectrum) . 

Using Eq. (|25[) and assuming the pair occupation n, to 
be uncorrelated, we have 

^ = 4]>>-A) 2 < . (A.l) 

i 

where ct 2 . is the variance of the pair occupation m. To 
demonstrate the validity of Eq. (|A. 1|) . we compare in 
Fig. [5] the exact intrinsic variance of n (solid circles) with 
the r.h.s. of Eq. (|A.1|) (open circles). These quantities 
can be calculated directly in CSMC. The results shown 
are for a uniform single-particle spectrum with a level 
spacing of 1 MeV at half filling and a constant pairing 
interaction of Vij = 0.3 MeV using a time step of v = 0.1. 
We have checked that Eq. (IA.1I) remains a good approxi- 
mation for the intrinsic error in more general cases, e.g., 
away from half filling and for a non-uniform spectrum. 

Since rii can take only two values or 1, the variance 
a\. of a given pair occupation is completely determined 
by its average 

al. =n;(l-ni). (A.2) 
Next, we use the BCS wave function to estimate 



where m and Vi are the usual BCS amplitudes. Combin- 
ing Eqs. (|A.ip . (IA.2I) and (|A.3|) . we have 

^=4^( £i -A) 2 r J^ 

= 25>( £l -A)^-^, (A.4) 
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FIG. 8: The intrinsic error cr m as a function of SI for a uni- 
form single-particle spectrum with level spacing of 1 MeV 
and a constant pairing interaction Vij = 0.3 MeV. The exact 
intrinsic error (solid circles) calculated directly from CSMC 
(using a time step of v = 0.1) is compared with the estimate 
of Eq. (|A.1|) with the exact a\ . calculated in CSMC (open cir- 
cles). The solid squares describe the BCS estimate Eq. (|A.4|) 
while the open squares correspond to Eq. ()26|) . The various 
lines are fits describing a scaling of fi 3 / 2 



where Aj are the level-dependent BCS pairing gaps. The 
BCS estimate (|A.4[) for a- m is shown in Fig. [5] by solid 
squares. 

The BCS expression (|A.4[) for the intrinsic variance can 
be further simplified for a uniform single-particle spec- 
trum (in which case the gap A is level-independent) with 
a bandwidth E c 3> A. In this limit most single-particle 
levels satisfy |e» — A| 3> A and the BCS amplitudes can 
be replaced by their non-interacting values. This leads 
to the simple expression in Eq. (j2"6"j) . 

In Fig.|S]we also compare the BCS estimate (|A.4I) (solid 
squares) with its simplified version (|26l) (open squares). 
The latter slightly overestimates the BCS expression; 
both underestimate the exact intrinsic variance but pro- 
vide a reasonable estimate within a factor of ~ 2. All 
cases in Fig.[8]scale as f2 3 / 2 (the respective fits are shown 
by lines). 

In the remaining part of this Appendix we show that 
Eq. (jA.ip (with A being the chemical potential) is a good 
approximation, i.e., that the covariance contribution to 
of n in (IA.1I) is small. 

Eq. (|25[) holds for an arbitrary parameter £ replacing 
A. Since the number of particles N is conserved, we have 
in general 



<4 = 4$>-C) 2 <+A(C) 



(A.5) 



where 



A (0 = 4 ^(e* - Ofo ~ C)cov(n i ,n i ) , (A.' 



and cov(ni, rij) = TiiUj — riifij is the covariance of rii and 
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FIG. 9: The variation of the fractional error A/af n (see Eq. 
()A.6|I ) as a function of £ for a uniform single-particle spectrum 
(level spacing of 1 MeV) with f2 = 20 at half filling and a 
constant pairing interaction Vij = 0.3 MeV. 



The covariance contribution to (|A.5|) vanishes for val- 
ues of C for which A(() — 0, a quadratic equation in 
£. Using particle number conservation 2 ^ . rii = N = 



2 52 . rij , we have 



Eq. (|A.7|) implies that 

2 



y^cov(nj,nj) 



(A.7) 



(A.8) 



Using Eq. (|A.8|) . we can express the coefficients of C 2 
and £ in the quadratic equation A(£) — in terms of the 
variances alone. We then find the following two solutions 



c± = Co ± 



1 + 



y~l£i£jCov(ni,nj) 

&3 



C 2 I> 



2 



(A.9) 



where 



Co 



(A.10) 



is the midpoint between the two solutions. 

As an example we show in Fig. [3] the quantity A/af n 
as a function of £ for a uniform single-particle spectrum. 
In general the zeros £± of A are not known without per- 
forming a full CSMC calculation. However, we observe 
that \A/af B \ is rather small in the region between and 
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£+ as compared with its typical value outside this region. 
Thus taking ( m ( in Eq. (|A.5[) and ignoring A leads to 
a good approximation to of n . 

The sums on the r.h.s. of Eq. (|A.10|) are dominated by 
those levels i for which hi is close to 1/2, i.e., by levels in 



the vicinity of the chemical potential A. Thus we expect 
£o to be in proximity to the chemical potential. We can 
also estimate Co directly from Eq. (|A.10|) using the BCS 
expressions for hi in Eq (|A.3[) . 
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